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ABSTRACT 


In the numerical modeling of the crack propagation in dynamic fracture using stationary elements, 
a discrete and sudden release of node at the crack tip creates spurious oscillation in the kinetic and 
strain energy values. In order to reduce the oscillation, a moving node element was utilized. This 
element can model a continuous crack tip movement more closely. The moving node element is 
compatible with surrounding regular isoparametric elements and no remeshing is required during the 
crack propagation. In addition, two different central difference schemes were compared, and their 
results were almost the same. 
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I. INTRODUCTION 


Inherent flaws in engineering materials can lead to a 
catastrophic failure due to unstable crack propagation. By 
eliminating the conditions and/or manufacturing defects 
(i.e., overloading, fabrication) which can lead to crack 
initiation, catastrophic failure can be prevented. In many 
structural components absolute prevention can not always be 
guaranteed. For such structures, catastrophic failure can be 
minimized by a crack arrest system. To ensure the proper 
incorporation of a crack arrest system, the effects of a 
propagating crack must be understood. These effects can be 
understood through the study of dynamic fracture mechanics. 

Dynamic fracture mechanics can be applied to any problem 
involving a body containing a crack in which inertia forces 
play an important role. In practice two kinds of dynamic 
fracture problems have received the most attention [Ref. 1]. 
These are: 

1. bodies with stationary cracks that are subjected to a 
rapidly varying load, and 

2. bodies under fixed or slow varying loading that 
contain rapidly moving cracks. 
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Dynamic crack propagation can be divided into two 
categories: crack initiation and crack propagation. A third 
category sometimes used is crack arrest. There are different 
opinions concerning crack arrest. Kanninen [Ref. 2] asserted 
that crack arrest is not a separate category of dynamic crack 
propagation, but is the termination phase of crack 
propagation. 

Dynamic fracture mechanics problems have focused mainly on 
bodies that contain a rapidly moving crack. As such, several 
numerical dynamic fracture analyses for the opening mode crack 
propagation have been studied. Kobayashi, et al. [Ref. 3] 
analyzed two fracturing Homalite-100 plates by dynamic finite 
element and dynamic photoelastic analyses. These analyses used 
the process of discrete crack-tip advances. The restraining 
nodal force was suddenly released when the crack-tip reached 
the next adjacent node. As indicated by Malluck and King [Ref. 
4], the above procedure has inherent problems: the sudden 
release of a node can induce unwanted high-frequency of 
motions; and the crack tip location within the nodal spacing 
can not be determined. To reduce this problem Malluck and King 
incorporated a mechanism for energy release in their finite 
element analysis. The nodal reaction of the crack tip was 
gradually reduced as the crack tip propagated (in the 
continuum view) from one node to the next node. 

Kobayashi, et al. (Ref. 5] were aware of the considerable 
oscillations in the calculated dynamic energy release rate. 
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They believed the oscillation was due to the sudden release of 
the crack tip (finite element nodes). To reduce these effects 
a nodal force release mechanism was incorporated to depict a 
more gradual transit of the crack tip between adjacent nodes. 
However, the assumed nodal release force is a somehow 
arbitrary choice. 

As indicated in [Ref. 6], computational procedures for 
crack propagation problems can be grouped into two types: 
stationary mesh procedure and moving mesh procedure. The above 
analyses may be categorized as a stationary mesh procedure. To 
predict the continuing propagation of a crack in a discrete 
model closely, moving mesh procedures were introduced. 
Nishioka and Atluri [Ref. 7] analyzed the propagating crac“L 
using a moving singular element in which the geometry of the 
element was altered as the crack propagated; thereby, 
requiring remeshing of the elements after an elapse period of 
time. If the element is distorted severely the accuracy may be 
degraded as indicated in [Ref. 8]; therefore, it was critical 
that remeshing occurred at the appropriate time. Kwon and Akin 
[Ref. 9] developed a procedure for modeling the crack 
propagation problem using a node moving along the edge of the 
element. In this procedure the element is not distorted; 
therefore, remeshing is not required. A modification of the 
latter procedure was incorporated into this study. 

The main objective of this study is to examine different 
numerical modeling techniques for studying dynamic cracJi 
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propagation. The main emphasis is placed on reducing the 
spurious oscillation in the calculated energy terms which 
resulted from the stationary node procedure. Figure la. An 
eight noded regular element (Figure lb) was used in the 
procedure. To reduce the spurious oscillation a moving node 
procedure (Figure 2a) using a moving node (Figure 2b) was 
incorporated into the finite element model. 

The study focused on the analysis of a rapidly propagating 
crack of opening mode in a linearly elastic isotropic body. 
Inertia was formulated into a diagonal mass matrix and the 
numerical time integration was accomplished using the two 
different central difference method. Crack propagation was 
simulated by the sequential release (at prescribe time 
intervals) of the nodes along the edge of the finite element 
model. Both stationary and moving node techniques were used to 
model the crack tip movement. Crack tip stress singularity was 
not represented in the model. From these method, crack opening 
displacement, work, strain energy, and kinetic energy were 
calculated and a comparative analysis was conducted from the 
results. 
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Figure 1. (a) Modeling of a propagating crack using a 

stationary node element, (b) An eight noded regular 
element. 






Figure 2. (a) Modeling of a crack propagating using a 

moving node element, (b) An eight noded moving node 
element. 







II. MATHEMATICAL DERIVATION 


A. GOVERNING EQUATION 

In an elastodynamic problem of a continuous isotropic 
body, certain governing equations must be satisfied at all 
interior points of the body. For two-dimensional problems, the 
equations of motion are: 


X ^ 

dx dy ^ dc 


^F-p 


d^u 

d^c 


= 0 


dv’ „ ^ r\ 

_-D-=0 

dx dy ^ dt y ^ dc^ 


( 1 ) 


where <jy, and r^y are the normal stresses in the x and y 
direction and the shear stress, respectively. and Fy are 
the body forces in the x and y direction, respectively, p is 
the density of the material, /x is the damping coefficient for 
the viscous damping, u and v are the displacements in the x 
and y direction, respectively. And t denotes the time. 

The constitutive equation (stress-strain relations) is 

{a} = [£?]{€} (2) 

where {a} represents the stress vector and {e} the strain 
vector, as shown below: 


7 





(3) 


{O} = [O, Oy Vj’’ 

{€} = [€^ 6 , 


and [D] is a symmetric 3x3 material property matrix. For an 
isotropic material, [D] is 


[D] 


1 V 


E 

(1-V-) 


V 1 
0 0 


0 

0 


1-v 


2 


(4) 


for the plane stress condition and 


[Z?]=. 


£:(i-v) 


(1+v){l-2v) 


1 

V 


1-v 

0 


1-v 

1 


0 


0 

0 

l-2v 


2 (1-v) 


(5) 


for the plane strain condition. 

The kinematic equation (strain-displacement relation) is 


•y 

tYxyJ 




du 
dx 
dv 
dy 
du ^ dv\ 
dy 


( 6 ) 
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B. BOUNDARY CONDITIONS 


In addition to satisfying the 90 verning equations, the 
prescribed conditions on stress and/or displacement must be 
satisfied on the boundary surface of the body. These 
conditions are classified as below. The essential boundary 
conditions are 


u=u 


v-v 


(7) 


and the natural boundary conditions are 

xy^x"*^ ^y^y 

where and <py are the traction in the x and y direction, 
respectively. And n^^ and n^ are the components of the outward 
directed unit normal vector on the boundary surface in the x 
and y direction, respectively. Boundary value problems, in 
general, consist of a combination of the two types described 
above. 

To transform the differential equations to matrix 
equations, apply the Galerkin method of weighted residual to 
the equations of motion (1): 
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(9) 



dx dy dc 



jo), dQ =0 



dy ^ dy 



(j^dQ =0 


where and u >2 are test (weighing) functions and is the 
domain of the given problem. Because Galerkin's method yields 
a non-symmetric matrix equation, apply the weak formulation to 
equation (9) to get a symmetric matrix. The results yield a 
symmetric equation because the governing equation is self- 
adjoint. Integrating by parts the stress terms of equation (9) 
yields: 




where T is the boundary surface of the domain n. Rewriting 
equation (lO) into matrix formr 




( 11 ) 
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where the first term becomes the stiffness matrix, the second 
term becomes mass matrix, the third term becomes the damping 
matrix, the fourth term becomes the load vector, and the fifth 
term becomes the body force vector. Substituti.ng the 
constitutive equation (2) into (11) gives 




0 

(i)o 





dx 

0 

8(i)j 

dy 

[V]' 




0 ' 

d^u 

dt^ 


0), 

0 ' 

rau) 

dc 

fa' 

0 

d(i). 



-+p 

0 


d^v 


0 

U)2, 

dv' 


dy 

dx 


Yxy, 



dt^ 




[dc 


>) do) = 


( 12 ) 


Substituting the kinematic equation (6) into (12) gives 


f 6<*), 



1 

du 







dx 

i 

d^u' 

1 . 

du 

dy 


dv 


<i)i 0 

dt‘ 


l"x 0 

dt 

0(*>2 d<j}2 

CO]. 

dy 

♦p 

3 

o 

d^v 


3 

o 

dv 

dy dx \ 


du ^ dv 


dt^J 


dc 

1 

dy dx 







(13) 


The matrices and the load vector are computed on the element 
level by discretization of the domain £2 and the boundary 
surface F; 


Q=EQ^ 

r=2:r* 


(14) 
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where Qg and Tg are the element domain and element boundary 
surface, respectively. 


C. MASS MATRIX 

The mass matrix comes from 


f 

“1 

0 ' 

Jq« 

0 

“2 


d^u 


dt^ 

d^v 


at" 


dQ 


(15) 


In this study eight noded isoparametric elements were used for 
the finite element mesh. The finite element derivation for 
the eight noded element is shown in detail below. The element 
is shown in Figure lb. 

From variable interpolation, 

u=Hj_u^^H 2 U2 U3 +H5 U5 +//s Ug U7 +//g Ug 

(16) 

v=H^ +H2 V2 +H3 V3 Vg +//5 V5 Vg V, Vg 


where u^ and v^ are the displacements in the x and y direction 
of the respective node and are the shape functions at the 
respective node. 

For the regular eight noded isoparametric element 
[Ref. 10] they are 
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//,=— (1-s) (1-c) 

’ 4 

H, = — (1-s) (1-C) (s-c-1) 

^ 4 

//, = — ( 1 -s) (1-c) (1-c) (s+c-1) 

’ 4 

H.=— (l-s) (1-C) ( c-s-1) 

(17) 

^5 = ^ (1-C) (1-S‘) 

W,= — (l-s) (1-c^) 

® 2 

//, = — (1-C) (1-s^) 

2 

//. = — (1-S) (1-cO 
* 2 


and for the eight noded isoparametric moving node element 
[Ref. 9] they are 


= ^ (1-S) (1-t) 

//2 = 4 (l^S) (1-t) 

4 Z Z 

H, = ^a*s) (l^t)-|w,--|H, 
H, = ^ (1-s) (l-C) 

//, =- - -(1-S‘) (1-t) 

2 (1-So) 

//6 = -| (1+S) (l-t^) 
fi, = ^ (1 + t) (l-s^) 

/fa = -| (1-s) (1-t^) 


(18) 


Since the shape functions are independent the time, the 
acceleration terms are 
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u. 

<y u 

dt^ 0 W2 0 w, 0 0 H5 0 0 H, 0 //, 0 I : 

d^v 0 0 Hx 0 H, 0 0 //, 0 0 //, 0 H, I : 

dC^ u. 


(19) 



Since the shape functions are expressed in the natural 
coordinates, it is necessary to carry out the integration in 
equation (21) in the natural coordinate, using the 
relationship 

dQ=dxdy=det [J] dsdt (22) 
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where [J] is Che Jacobian macrix given as 

dx By 
Bs Bs 
Bx By 
Be Be 


SubstiCuting equation (22) into (21) yields 

u, 

V; 

U5 
V; 

dec [ J] dsde ' 

Since the shape functions are an explicit function of s 
and t, a numerical method is used to evaluate equation (24). 
The two-point Gaussian quadrature formula is used to carry out 
Che integration. The result is a symmetric consistent mass 
matrix. 

For this study the consistent mass matrix was diagonalized 
by following procedure [Ref. 11]: 

a. Calculate the diagonal coefficients of the consistent 
mass maCrix. 



0 


0 


0 

0 


0 

.. 

0 



0 


0 

. . 

0 

0 


0 

.. 

0 



0 


0 

.. 

0 

0 


0 

.. 

0 



0 


0 . . 

.. 

0 

0 


0 

H,H, .. 

0 

^9^8. 




(23) 
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b. Determine the total mass of the element, m. 

c. Add the diagonal coefficients, associated with 

translation only to obtain a number, e. 

d. Multiply the diagonal coefficients, by m/e. 

e. Set all the off-diagonal terms to zeros. 


D. DAMPING MATRIX 

The damping matrix is given by 



0 


du 


dt 

dv 


dt 


dQ 


(25) 


Since the shape functions are independent of time. 



du 








dt 


0 

H2 0 


0 


0 //g 0 Hg 0 H, 0 

dv 

[at] 


0 

0 H2 

0 

^3 

0 

0 

0 

0 


//g 0 

0 H, 


(26) 


Substituting equations (20), (22) and (26) into (25) and 

applying the two-point Gaussian formula gives a symmetric 
consistent damping matrix. 

E. STIFFNESS MATRIX 

The stiffness matrix is given as 
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d(j>. 

dx 


0 



5<i), 

dy 


8(0, 


dy 

3 ( 0 - 


dx 


[D] 


du 

dx 

dv 

dy 


}dQ 


du ^ dv\ 

dy 8 xJ 


(27) 


Rewriting the kinematic equation in terms of nodal 


displacements gives 


du 

dx 


I ^ \ 

dy 

du ^ ^ 
dy dx 



0 


dy 


ib 

dx 


' dH^ 

0 

dH^ 

0 . . 

r \ 

dH^ 


8W3 

0 

u.; 

dx 

dx 

. . u 

dx 

u 

dx 

V. 

0 

dH^ 

dy 

0 

dH^ 

dy 

. . 0 

dH^ 

~dy 

0 

dH, 

dy 


‘ 

3H^ 

dH^ 

dH, 

dH^ 

dfu, 

dH^ 

8H3 

dh’g 



~dy 

dx 

~d^ 

dx 


dx 

dy 

dx 




(28) 


where 


dx “ ds 


A 

dt 

dH, 




dt 


(29) 


and A^j are the components of the inverse of the Jacobian 
matrix. The matrix in ecjuation (20) is called the [B] matrix 
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and the vector is called the displacement vector (d} . From 
Galerkin's method on the first matrix in equation (27) gives 


3(i), 

dx 

0 


0 

dy 


8(0. 

dy 

8 ( 0 ^ 

dx 


[B] ■ 


(30) 


Substituting equations (29) and (30) into (28) gives 


f [BldQid] 

J Q- 


(31) 


The stiffness matrix is given as 

[icn=f [B]^lD}[B]dQ (32) 

J Q« 

Substituting equation (22) into (32) gives 

[iC^] [fl] ^[£?] [Bidet [^] dsdc (33) 

The two-point Gaussian quadrature formula is used to carry out 
the integration, 

P. BODY FORCE VECTOR 

The body force vector is given as 
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Substituting equations (20) and (22) into (34) gives 


iH, 0 Q 

-iJ-i 0 0 



0 He 

H, 0 


0 

^9 



det [ J] dsdc 


(35) 


The two-point Gaussian quadrature formula is used to carry out 
the integration. 


G. LOAD VECTOR 

The load vector is given by 



(36) 


where F is the boundary surface of the domain 0. If there is 
no traction 



(37) 


If there is traction, substituting equation (20) into (36) 
gives 



- [jr 0 JT; 0 r 

{p} =r _ _ _i ih) ds 

J4-7-3 0 0 H, 0 H.. 2 


where {P} is Che external load vector along edge 4-7-3 (Figure 
3) , h is the element thickness and overbars indicate the shape 
functions are evaluated at t=l (the natural coordinate). 
Integrating equation (38) yields 


{P}=. }h 


1 i) 


Fraction of the total force is shown in Figure 3. 


Ty- - 2/3 

Oy- 0,- 1/6 


Figure 3 . Fraction of total 
force on the element 
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From the procedural derivation above, the equation of 
motion is finally rewritten as 

[M] {d}>[C] {d} = {;?} (40) 

where {r}={p}+{f}. This equation applies to any dynamic 
problem in which the finite element method is used. 

H. TRANSIENT ANALYSIS 

In dynamic problems the displacements, velocities, 
accelerations, strains, stresses, and loads are time 
dependent. Therefore, a time integration scheme is required. 
The equation of motion at time t is 

[M] {d]^={R}^ (41) 

This study examined two different forms of the central 
difference method: 

The first method is as summarized below [Ref. 12]; 

A. Initial Calculations 

I. Solve for [M], [C], and [K] 

2. {d}'' = [M]-i[{R}0 - [C]{d}0 - [K]{d}°]. 

where {d}° and {d}° are the initial conditions. 

3. {d}-^*^ = {d}0 - (At){d}0 + a3{d}0. 

4. [M] = ao[M] + ai[C]. 
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B. For each time step 

1. {R} = {R}’' - ( [K] - a2[M]){u}'' - (ao[M] - a^ [C] ) 

2. = [M] 

3. {d}*^ = ao[{d}t"^‘^ - 2{d}'^ . {d}"-^*^] . 

4. {d}*^ = - {d}'^-^'^] 

where ag = 1/At^; a^ = l/2At; a 2 = 2ao; a 3 = l/a^ 

The second form "summed form" [Ref. 13] is as follows: 
A. For each time step 

1. {d}^'^-"^ = [MI'Mr^'^*^ - [K] {d}^'^"’M . 

2. ^ |^|At + l/2 ^ + 

3. {d}^‘'-"2 , {d}^t + l ^ ^tAt-H3/2|^}At + 3/2^ 


The explicit central difference scheme is conditionally 
stable. The time step At must be controlled from the smallest 
period of the system. As indicated in [Ref. 14], a safe 
limiting value is given by 


At=0.45L" 


p (l-»-v) (l-2v) 
E{l-v) 


1/2 


(42) 


where L" is the smallest distance between adjacent nodes of 
any element used .'>rid E and p are the elastic modulus and 
poisson's ratio, respectively. This safe limiting value was 
used as a criteria to ensure stability. 
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I. DYNAMIC FRACTURE ANALYSIS 

From elasticity the stresses near a crack tip for the 
opening mode are: 



-^cos—(l -sin — sin —j 

2WI 2 2 I 


v'2n 

K 


42 


K 6/, . 0 30 

— cos— 1 +sin— sin — 

Tir 2\ 2 2 


K I ■ & 6 30 

—-—sin —cos —cos — 

v/ItT? ',222 


(43) 


where and Cy are the tensile stresses in the x and y 
direction, respectively, is the shear stress, r is the 
radial distance from the crack tip, 0 is the angle from the 
crack plane, and K is the stress intensity factor. The stress 
intensity factor provides a convenient parameter of the stress 
distribution around a flaw. If K is equal to or exceeds a 
critical value (fracture toughness, K^.) the crack becomes 
unstable. is a material property. [Ref. 15] 

In the dynamic generalization of linear elastic fracture 
mechanics (LEFM) , has two parts. For the initiation of 
crack growth 

K{a,o,t) = K^{o) (44) 


where a represents the crack length, a is the applied stress, 
t is time, and a is the loading rate. Similarly, for a 
propagation crack 
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Kia, a. t) = K^iA) 


( 45 ) 


where a denotes the crack speed. [Ref. 2] 

An alternative criterion for determining the motion of a 
crack is, 

G(a,o, t) = i?(i) (46) 

where R is the energy dissipation rate required for crack 
growth and G is the dynamic energy release rate given by, 

G = i/ dW_ dU _ dT] 
b\da da da) 

. 1 I dU_ dT) 
bi\dt dt dtj 


where U is the strain energy, T is the kinetic energy, W is 
the work done on the structure by external loads, a is the 
crack length, and b is the plate thickness at the crack tip. 

The dynamic energy release rate and the dynamic stress 
intensity factor are connected through Freund's formula. 


Ko 


\ EG 

U-v" 


AU) 



(48) 


where G is the dynamic energy release rate, E is Young's 
modulus, p is Poisson's ratio, and A is a function that is 
dependent on crack speed. A is given as 
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where Cj^ and C 2 are the longitudinal and shear wave speeds. 
[Ref. 2] 

The wave speeds are given by [Ref. 16], 


c. 


P 


(50) 


where X and pL are Lame's constants, and 


25 




III. RESULTS AMD DISCUSSION 


A. PROBLEM STATEMENT 

The main objective is how to best approximate the behavior 
of a continuous propagating crack using a discrete finite 
element model. The model was composed of a uniform mesh of 
eight noded isoparametric elements. In the first model the 
nodes were stationary. To approximate the movement of the 
crack, the nodal restraints were suddenly released at given 
time intervals (time required for the crack to travel a nodal 
spacing in the continuous movement) . As will be shown, this 
resulted in spurious oscillation in the energy terms. To 
reduce the oscillation due to discrete jumps of the crack 
movement, as simulated by the model, a moving node element was 
employed in the second model. 

Each of the analysis was based upon the plane-strain 
conditions of a centered cracked plate of homogeneous 
isotropic elastic material. The plate deformed under the 
action of an applied time-independent (constant) traction. The 
applied load was normal to the crack plane, mode I 
configuration. The mass matrix was diagonalized and the 
numerical time integration was accomplished by two different 
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forms of the central difference method. A comparative study 
was conducted of the two forms. 


B. VERIFICATION EXAMPLES 

The developed finite element analysis code was verified by 
analyzing a problem with a known solution. The problem 
considered is governed by the wave equation. A long uniform 
rod with one end fixed and the other end free was -xamined. 
The free end was subjected to an applied axial force which was 
released at time zero. Applying Newton's second law and for a 
constant area multiply by Young's modulus (AE) the equation of 
motion simplified to the wave equation 


dC^ dx^ 


(51) 


where c^ = E/p and u is the axial displacement of an element 
of the rod. By the method of separation of variables, the 
analytical solution is expressed as 


uix, t) = 


8FL 


TT 


(- 1 )" 
(2n+l)2 


sin (2n^l)7ix ^P3J2n^l)7rct 


(52) 


2L 


2L 


As shown in Figure 4, the numerical result was reasonably 
accurate. To verify the accuracy of the second model (moving 
node), it was compared to the stationary model when node five 
was at the center edge of the element. For this case the 
results should be and were the same. 
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Figure 4 . Displacement of end of beam 
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As another comparative check of the accuracy of the 
developed methods, the results of the crack openina 
displacement (COD) was compared against the analytical 
solution as derived by Broberg [Ref. 17] . To obtain the 
mathematical description of a two-dimensional crack, Broberg 
made the foll.,wing assumption [Ref. 17] : 

1. A two-dimensional crack can propagate in one plane only. 

2. The material is homogeneous and isotropic with regard to 
fracturing characteristics and stress-strain relations. 

3. The material obeys Hooke's law and is perfectly brittle. 

4. The surface energy is zero. 

5. The crack propagates with a maximum velocity from the 
start. 

Broberg's analytical solution applies to an infinite medium 
with the above characteristics. These assumptions were 
incorporated into this study. The comparison between 
analytical and numerical solutions are shown later. 

C. NUMERICAL TIME INTEGRATION 

An explicit procedure based on the central difference 
method was incorporated into this study. The two forms of the 
method are as outlined in Chapter II. The latter is the 
"summed form" [Ref. 13]. From a computational consideration, 
the "summed form" has an advantage. For a damped structure, 
the "summed form" only requires the mass matrix be 
diagonalized to simplify the computation; whereas, the former 
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requires the diagonalization of the mass and damping matrix 
for simplification. To date no supporting argument or test 
cases exist which justifies the use of a diagonal damping 
matrix. 

During the study it was noted that the first central 
difference method was very sensitive to computational round 
off errors. As the time increment decreased, t.he central 
difference did not converge while running the model in single 
precision. The model was ran at double precision (minimizing 
round off error) to achieve convergence. On the other hand, 
the "summed form" converged during single precision runs. 
Comparison of the convergent results revealed no difference 
between the two central difference techniques. The "summed 
form" simplifies the computation, and as a result, it saves 
time and money. 

An important consideration when using the central 
difference method is the time step. Because the central 
difference method is conditionally stable, the time step must 
be less than the critical time step to ensure stability. The 
critical time step is defined as the time required for an 
acoustic wave to transverse the smallest element of the system 
[Ref. 11]. A time step less than the critical time step 
guarantees the stability of the solution. Stability meaning 
that the displacement does not grow without limit. As 
indicated in [Ref. 11], the practical limit on At is 
approximately 25% less than 2/ci}^^. As seen in this study. 
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decreasing the time step size considerably below the 
approximate formula for At^^., equation (42j, resulted in no 
change. 

In the second analysis, the distance the moving node can 
travel is limited by /i. H is the fractional distance of the 
element length. It represents the minimum distance the moving 
node should be from the corner node. Figure 5. Because 
convergence had occurred, decreasing At for a constant beta 
resulted in no change in the solutions. The minimum distance 
from the corner node, 0, is dependent on At for relatively 
large time step size. As shown in Figure 6, decreas.*.ng 0 for 
a constant At introduced higher frequencies with similar 
magnitudes. If the moving node propagated within the minimum 
allowable distance, the system became unstable. Exceptionally 
high frequencies and amplitudes resulted, Figure 7. 

As shown in the Table 1, further decreasing the time steps 
resulted in a constant 0. 


Table 1. Fractional minimum distance of 
the moving mid-node from the end node 


NSTEP 

10 

25 

50 

100 

beta {0) 

1/5 

1/10 

1/10 

1/10 


*Nstep is the number of time steps for 
the crack to advance one half the element 
length. 
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to element length 


D. STUDY OK MASS AND STIFFNESS MATRICES 

As indicated in [Ref. 11], test cases have shown that the 
diagonal matrix, using the procedure outlined in Chapter II, 
gives greater accuracy than a consistent mass matrix. A 
diagonal mass matrix also requires less storage space than a 
consistent mass matrix and simplifies the computation of the 
central difference scheme. 

For the regular eight noded isoparametric element, the 
concentrated mass at the corner nodes are one-eight that of 
the mid-nodes. For the moving nodal case, the diagonalized 
mass matrix from the moving node element caused unstable 
(oscillatory) conditions. Because of the close proximity of 
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moving node five to node one, the diagonal mass matrix 
procedure resulted in the majority of the mass concentrated at 
nodes five and one. As node five moved across the edge of the 
element toward the center, the nodal point mass approached 
that of the regular element. Due to the unstable conditions 
resulting from the diagonal mass procedure for the moving node 
element, the diagonal mass matrix was calculated with the 
moving node at the center of the element boundary edge. The 
computed diagonal mass matrix remained constant throughout the 
analysis. 

For the eight noded regular element the element nodal 
stiffness remained constant. On the other hand, for the eight 
noded moving node element the nodal stiffness of nodes one, 
two, and five varied as node five moved along the edge of the 
element. Node five cared the higher stiffness. The remaining 
nodes (three, four, six, seven, eight) were equivalent to the 
regular element. 

E. COMPARISONS OF ENERGY TERMS 

In order to make a direct comparison between the use of a 
stationary element and a transition element, the generated 
results were compared with Broberg's problem. The finite 
element breakdown for Broberg's problem was as depicted by 
Kobayashi et al. [Ref. 5] and as shown in Figure 8. 

Figure 9 shows the crack opening displacement during crack 
propagation. The two finite element results are compared to 
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trhe analytical solution as solved by Broberg. Except for some 
positions, the results were relatively similar. No comparative 
conclusion can be drawn from the results. Although their close 
comparison to the analytical solution indicates the 
reasonableness of the model. 

Figure 10 shows the crack opening displacement at x = 0. 
From this result, it is seen that the use of the moving node 
element approximates a propagating crack more closely than the 
stationary element. 

Figure 11 shows the work on the plate during crack 
propagation. There is a similar trend between the stationary 
element and the moving node element with the exception that 
there is a slight increase in the work for the moving node 
element. 

Figure 12 shows the calculated strain energy during crack 
growth. The stationary element produces spurious oscillation. 
The oscillation is caused by the instant release of the crack 
tip during the process of discrete crack tip advances [Ref 5] . 
As shown in Figure 12a, there was a rapid build up in the 
strain energy until the node was released. Upon nodal release 
there was a rapxd reduction in the strain energy followed by 
a rapid increase until the next nodal release. The oscillation 
amplitude of the strain energy increased during crack 
propagation. Employing a moving node element made the crack 
tip movement more continuous. As a result, the oscillation 
amplitude of the strain energy was much reduced. Figure 12b. 
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Figure 13 shows the calculated kinetic energy during crack 
growth. The delayed nodal release, caused by the stationary 
element, resulted in a rapid decrease in the kinetic energy 
prior to the release of the node, then an rapid increase 
followed by a rapid reduction until the next nodal release. 
This trend was a reversal of that seen for strain energy. The 
oscillation amplitude of the kinetic energy increased during 
crack growth as it did for the strain energy. Employing a 
moving node element reduced the oscillation amplitude because 
the moving node better represents a moving crack tip. 

Although not necessary since the crack had already passed, 
if variable interpolation was used to move the middle node 
(node five) back to the center position, there was an increase 
in the oscillation of the results. The increase in oscillation 
is probably due to the abrupt change in the nodal stiffness 
when bringing node five back to its center location. As the 
crack propagated, there was a gradual shift in the stiffness 
of node five. 

In an elastic medium there are only two types of waves 
that propagate through a solid that is unbounded [Ref. 16]. 
These are dilatational and distortional waves. If the solid 
has a free surface, Rayleigh surface waves may also propagate. 
Rapid reduction in stress occurred upon nodal release. Due to 
this rapid reduction, unloading waves were generated at the 
crack tip. As the wave propagated through the medium, strain 
energy decreased and kinetic energy increased. 
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If Che body is bounded, elastic waves reflect and refract 
off Che boundaries back toward the crack tip as loading waves. 
The loading wave expands the body; thereby, amplifying the 
strain caused by the tension load. This result of an increase 
in strain energy is seen in Figure 13a. 
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Figure 6. Effect of varying beta on strain energy for 
nstepl = 50: (a) beta = 1/20, (b) beta = I/IO, (c) beta 
1/5. 
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Figure 8. Finite element mesh for Broberg' 
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Figure 9. Crack Opening Displacement for (a) stationary 
element and (b) moving node element 
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Figure 10. Crack opening displacement at x = 0 
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Figure 11. Work vs Crack Length (a) stationary element, 
(b) moving node element (beta=l/5) 
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Figure 12, Strain energy vs Crack Length for (a) 
stationary element (b) moving node element {beta=l/5) 
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IV. CONCLUSIONS AND RECOMMENDATIONS 


This study focused on the finite element modeling of a 
rapidly propagating crack. Two methods were empxoyed in this 
study. The first procedure, adopting stationary elements, 
resulted in spurious oscillation of the energy terms. For this 
procedure the propagating crack was represented by discrete 
jumps at given time intervals. To more closely approximate the 
propagating crack, a moving node procedure was employed in the 
model. As shown in the previous chapter, the moving node 
procedure reduced the spurious oscillation of the energy 
terms. Unlike the moving mesh procedure, no remeshing is 
required for this procedure. 

From the comparison of the two procedures mentioned above, 
the following observations were noted: 

1. The two forms of the central difference methods yield 
almost the same results. As a result, the "summed form" was 
preferred due to its simplicity and efficiency in 
computational time. 

2. The range which the moving node, which represents the crack 
tip, can propagate along the edge of the element was a 
function of the time step size for relatively large time step 
sizes. However, for reduced time step sizes, the range was the 
same. 


45 






3. For the moving node procedure, diagonalizing the mass at 
each increment of movement resu.'ted in sporadic results. When 
node five is close to the corner nodes the majority of the 
concentrated element mass was located at node five and the 
respective corner node it was close to. 

This study did not account for the singularity in the 
crack tip. To represent the stress field near the crack tip 
more accurately, it is recommended that singular element, as 
derived in [Ref. 9], be incorporated into the model. 
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